# -*- coding: utf-8 -*-
"""
@authors: coude

This library is designed to smoothh the HAWC+ polarization data to the Herschel resolution
    
"""

# REMINDER: Python arrays are inverted [y,x] relative to IDL [x,y]

# =============================================================================
# Package dependencies
# =============================================================================

from astropy.io import fits
import math
import numpy as np
from reproject import reproject_exact
from astropy.convolution import convolve, Gaussian2DKernel
import matplotlib.pyplot as plt
from scipy import stats, optimize

# =============================================================================
# Function to reproject an HAWC+ data product to the Herschel pixel scale
# =============================================================================
def ReprojectCube(SourcePol, SourceProj, Smooth=False):
	# Smoothing if allowed
	if Smooth==True:	
        # Determining the required smoothing FWHM
		ConvertFWHMtoDev = (2.0*(2.0*math.log(2.0))**0.5)**-1.0
		ResolutionPol = SourcePol[0].header['BMAJ']
		print('Beam Size of polarization data')
		print(ResolutionPol*3600)
		ResolutionProj = SourceProj[0].header['BMAJ']
		print('Beam size of target pojection')
		print(ResolutionProj*3600)
		PixSizePol = SourcePol[0].header['CDELT2']
		print('Pixel size of the polarization data')
		print(PixSizePol*3600)
		PixSizeProj = SourceProj[0].header['CDELT2']
		print('Pixel size of the target projection')
		print(PixSizeProj*3600)
		Smoothsize = ConvertFWHMtoDev*(ResolutionProj**2.0-ResolutionPol**2.0)**0.5
		SmoothsizePix = Smoothsize/PixSizePol
		print('Target standard deviation for smoothing in arcseconds')
		print(Smoothsize*3600)
		print('Target standard deviation for smoothing in pixels')
		print(SmoothsizePix)
		SourcePol = SmoothStokes(SourcePol,SmoothsizePix)
	# Reprojecting all the extensions of the HAWC+ file
	NewCube = ReprojectExtensions(SourcePol, SourceProj[0].header)
	# Recalculating the polarization information
	FinalCube = recalculate_polarization(NewCube)
	return FinalCube

# =============================================================================
# Function to create a fits container with a single new header 
# =============================================================================
def ReprojectExtensions(hawcfile, RefHeader):
	# Creating multi-extension FITS file
	newfile = fits.HDUList()
	for i in range (0,16):
		# Reprojecting the file for a given header
		projfile, footprint = reproject_exact(hawcfile[i], RefHeader)
		# Creating the HDU exension
		listfile = fits.ImageHDU(data=projfile, header=RefHeader)
		# Appending current extension to the new fits file
		newfile.append(listfile)
	
	# Adding the original tables to the final file
	# These may need to be modified to account for the change in coordinates
	newfile.append(hawcfile[16])
	newfile.append(hawcfile[17])
	return newfile


# =============================================================================
# Function to calculate the polarization properties from Stokes I, Q, and U
# from a fits container opened from an HAWC+ fits file
# See HAWC+ white paper for structure (Gordon et al. 2018)
# https://ui.adsabs.harvard.edu/abs/2018arXiv181103100G/abstract
# !! The calculation below assumes zero covariance as information is missing !!
# =============================================================================
def recalculate_polarization(newfile):
    # Creating human-readable names for arrays
	StokesI = newfile[0].data # Stokes I
	dI = newfile[1].data # Error Stokes I
	StokesQ = newfile[2].data # Stokes Q
	dQ = newfile[3].data # Error Stokes Q
	StokesU = newfile[4].data # Stokes U
	dU = newfile[5].data # Error Stokes U

	# Biased polarized intensity 
	PI_biased = (StokesQ**2.0 + StokesU**2.0)**0.5
	# PI uncertainties
	dPI = ((StokesQ*dQ)**2.0 + (StokesU*dU)**2.0)**0.5/PI_biased 
	# Mapping where the intensity is smaller than the uncertainty
	pimask = np.where(PI_biased < dPI)
	dPI_debias = np.copy(dPI)
	dPI_debias[pimask] = 0.0 #np.nan
	# De-biased polarized intensity PI 	
	PolPI = (PI_biased**2.0 - dPI_debias**2.0)**0.5  
	PolPI[pimask] = 0.0

	# Polarization fraction P
	P_biased = 100.0*PI_biased/StokesI
	PolP = 100.0*PolPI/StokesI
    # Uncertainty of polarization fraction
	dP = np.absolute(P_biased*((dPI/PI_biased)**2.0 + (dI/StokesI)**2.0)**0.5)
	# Eliminating unphysical results (setting values to 0.0 or 100.0)
	P_biased[np.where(P_biased < 0.0)] = 0.0
	P_biased[np.where(P_biased > 100.0)] = 100.0
	PolP[np.where(PolP < 0.0)] = 0.0
	PolP[np.where(PolP > 100.0)] = 100.0
	dP[np.where(dP > 100.0)] = 100.0
    
    # Polarization angle O
	PolO = (0.5*180.0/math.pi)*np.arctan2(StokesU, StokesQ)
    # Uncertainties dO
	dO = (0.5*180.0/math.pi)*((StokesQ*dU)**2.0 + 
                              (StokesU*dQ)**2.0)**0.5/PI_biased**2.0
    # Polarization angle B
	PolB = PolO + 90.0
	PolB[np.where(PolB > 90.0)] = PolB[np.where(PolB > 90.0)] - 180.0

	# Creating new fits container with updated polarization in ICRS
	workfile = fits.HDUList()
	# Stokes I and dI
	workfile.append(fits.ImageHDU(data=StokesI, header=newfile[0].header))
	workfile.append(fits.ImageHDU(data=dI.data, header=newfile[1].header))
	# Stokes Q and dQ
	workfile.append(fits.ImageHDU(data=StokesQ.data, header=newfile[2].header))
	workfile.append(fits.ImageHDU(data=dQ.data, header=newfile[3].header))
	# Stokes U and dU
	workfile.append(fits.ImageHDU(data=StokesU.data, header=newfile[4].header))
	workfile.append(fits.ImageHDU(data=dU.data, header=newfile[5].header))
	# Image mask
	workfile.append(fits.ImageHDU(data=newfile[6].data, header=newfile[6].header))
	# Polarization fraction
	workfile.append(fits.ImageHDU(data=P_biased, header=newfile[7].header))
	workfile.append(fits.ImageHDU(data=PolP, header=newfile[8].header))
	workfile.append(fits.ImageHDU(data=dP, header=newfile[9].header))
	# Polarization angle
	workfile.append(fits.ImageHDU(data=PolO, header=newfile[10].header))
	workfile.append(fits.ImageHDU(data=PolB, header=newfile[11].header))
	workfile.append(fits.ImageHDU(data=dO, header=newfile[12].header))
	# Polarized intensity
	workfile.append(fits.ImageHDU(data=PI_biased, header=newfile[13].header))
	workfile.append(fits.ImageHDU(data=dPI, header=newfile[14].header))
	workfile.append(fits.ImageHDU(data=PolPI, header=newfile[15].header))
	# Adding catalog at the end
	workfile.append(newfile[16])
	workfile.append(newfile[17])

    # Returning FITS container with updated polarization values
	return workfile

# =============================================================================
# Function to smooth the Stokes parameters of a polarization cube
# =============================================================================
def SmoothStokes(PolCube,SmoothSize):
	# Smoothing the Stokes parameters, with SmoothSize in pixels
	gauss_kernel = Gaussian2DKernel(SmoothSize)
	# Stokes I and dI
	PolCube[0].data = convolve(PolCube[0].data, gauss_kernel)
	PolCube[1].data = convolve(PolCube[1].data, gauss_kernel)
	# Stokes I and dI
	PolCube[2].data = convolve(PolCube[2].data, gauss_kernel)
	PolCube[3].data = convolve(PolCube[3].data, gauss_kernel)
	# Stokes I and dI
	PolCube[4].data = convolve(PolCube[4].data, gauss_kernel)
	PolCube[5].data = convolve(PolCube[5].data, gauss_kernel)
	return PolCube

# =============================================================================
# Creating a mask to downsample the data for the vector catalogs
# =============================================================================
def DownsampleVectors(FitsFile, Step=1):
    # Creating the empty array for the final mask
    maskref = fits.PrimaryHDU(data=FitsFile[0].data, header=FitsFile[0].header)
    mask = maskref.copy()
    # Set all elements to NaN
    mask.data.fill(np.nan)
    # Size of the map
    MaskSize = np.shape(FitsFile[0].data) 
    for i in range (0,MaskSize[0],Step): # Latitude
            for j in range (0,MaskSize[1],Step): # Longitude
                mask.data[i,j] = 1.0 # Setting every pixel at a given number of steps to 1.0
    return mask

# =============================================================================
# Creating the polarization vector catalog
# =============================================================================

def MakePCats(Nref, Tref, Iref, Pref, CatMask=None, IdI=10.0, PdP=3.0, Pmax=30.0):
    Ndata = Nref.copy()
    Tdata = Tref.copy()
    Idata = Iref.copy()
    Pdata = Pref.copy()

    # Creating lists
    CatI = [] # Stokes I
    CatdI =[] # Uncertainty dI for Stokes I
    CatQ =[] # Stokes Q
    CatdQ =[] # Uncertainty dQ for Stokes Q
    CatU =[] # Stokes U
    CatdU =[] # Uncertainty dU for Stokes U
    CatP = [] # Polarization fraction P (debiased)
    CatdP =[] # Uncertainty dP for polarization fraction P
    CatO = [] # Polarization angle O
    CatB = [] # Rotated polarization angle B (+90 degrees)
    CatdO =[] # Uncertainty dO for polarization angles
    CatPI = [] # Polarized intensity PI (debiased)
    CatdPI =[] # Uncertainty dPI for polarized intensity PI
    CatN = [] # Column density
    CatT =[] # Dust temperature

    fluxconv = 1000.0*(3600.0*Idata[0].header['CDELT2'])**-2.0 # Conversion factor from Jy/pixel

    # Creating a new fits object for APLpy's FITSFigure method to recognize
    # for the polarization fraction P
    pref = fits.PrimaryHDU(data=Pdata[8].data, header=Pdata[0].header)
    pmap = pref.copy()
    # for the magnetic field angle B
    oref = fits.PrimaryHDU(data=Pdata[10].data, header=Pdata[0].header)
    omap = oref.copy()
    # Creating a mask to hide polarization vectors with low signal-to-noise ratios
    imask_01 = np.where(Pdata[0].data/Pdata[1].data < IdI) # Total intensity SNR threshold
    imask_02 = np.where(Pdata[0].data < 0) # Total intensity positive threshold
    pmask = np.where(Pdata[8].data/Pdata[9].data < PdP) # Polarization SNR threshold
    pmaxmask = np.where(Pdata[8].data > Pmax) # Polarization SNR threshold
    # Masking all the indices for  which the selection criteria failed
    pmap.data[imask_01] = np.nan
    pmap.data[imask_02] = np.nan
    pmap.data[pmask] = np.nan
    pmap.data[pmaxmask] = np.nan

    # If region mask is provided
    if CatMask != None:
        Maskdata = CatMask.copy() # Mask for the catalog
        cmask = np.where(Maskdata.data != 1.0) # Catalog mask
        pmap.data[cmask] = np.nan

    # Measuring the size of the new array
    shapeX = Ndata[0].data.shape[0]
    shapeY = Ndata[0].data.shape[1]

    for i in range (0,shapeX):
        for j in range (0,shapeY):
            if pmap.data[i,j] > 0.0: # Iterate only on pixels where the selection criteria are respected
                CatI.append(Pdata[0].data[i,j]* fluxconv)
                CatdI.append(Pdata[1].data[i,j]* fluxconv)
                CatQ.append(Pdata[2].data[i,j]* fluxconv)
                CatdQ.append(Pdata[3].data[i,j]* fluxconv)
                CatU.append(Pdata[4].data[i,j]* fluxconv)
                CatdU.append(Pdata[5].data[i,j]* fluxconv)
                CatP.append(pmap.data[i,j])
                CatdP.append(Pdata[9].data[i,j])
                CatO.append(Pdata[10].data[i,j])
                CatB.append(Pdata[11].data[i,j])
                CatdO.append(Pdata[12].data[i,j])
                CatPI.append(Pdata[15].data[i,j]* fluxconv)
                CatdPI.append(Pdata[14].data[i,j]* fluxconv)
                CatN.append(Ndata[0].data[i,j])
                CatT.append(Tdata[0].data[i,j])

    # Transforming lists into numpy arrays
    CatI = np.array(CatI)
    CatdI = np.array(CatdI)
    CatQ = np.array(CatQ)
    CatdQ = np.array(CatdQ)
    CatU = np.array(CatU)
    CatdU = np.array(CatdU)
    CatP = np.array(CatP)
    CatdP = np.array(CatdP)
    CatO = np.array(CatO)
    CatB = np.array(CatB)
    CatdO = np.array(CatdO)
    CatPI = np.array(CatPI)
    CatdPI = np.array(CatdPI)
    CatN = np.array(CatN)
    CatT = np.array(CatT)


    # Returning all the catalogs 
    return CatI, CatdI, CatQ, CatdQ, CatU, CatdU, CatP, CatdP, CatO, CatB,CatdO, CatPI, CatdPI, CatN, CatT

# =============================================================================
# Polarization fraction as a function of total intensity
# =============================================================================
def PvI(CatP, CatdP, CatI, Iscale=None, Pscale=None, showfit='false', weighted='true', errorbars='false', FigSize=[5.12,3.84], Source=None):  
    PvI_plot = plt.figure(figsize=(FigSize[0],FigSize[1])) # Creating the figure object
    plt.xlabel('Stokes Total Intensity $I$ (mJy arcsec$^{-2}$)')
    plt.ylabel("Polarization Fraction $P$ (%)")
    # Plots with and without errorbars 
    if errorbars=='true':
        markers, caps, bars = plt.errorbar(CatI, CatP, yerr=CatdP, fmt ='.', color='black', ecolor='gray', markersize=1, mfc='none',elinewidth=0.9)
        [bar.set_alpha(0.25) for bar in bars]
    else:
        plt.plot(CatI, CatP, '.', color='black', markersize=1, mfc='none')         
    # Switching the total intensity in log scale
    plt.xscale('log')
    axes = plt.gca() # Calling the axes object of the figure
    axes.xaxis.set_ticks_position('both') # Adding ticks to each side
    axes.yaxis.set_ticks_position('both') # Adding ticks to each side
    plt.tight_layout() # Using all available space in the plot window

    # Checking the automated limits for I and P
    xmin, xmax = plt.xlim()
    ymin, ymax = plt.ylim()
    # Setting mininum and maximum values of the I plot
    if Iscale==None:
        Iscale = np.empty(2)
        Iscale[0] = xmin
        Iscale[1] = xmax
    if Pscale==None:
        Pscale = np.empty(2)
        Pscale[0] = ymin
        Pscale[1] = ymax
    # Setting the minimum and maximum values
    plt.xlim(Iscale[0],Iscale[1])
    plt.ylim(Pscale[0],Pscale[1])


    if showfit=='true':
        # Create the array with ln(I), ln(P), and d(ln(P))
        # The numpy package uses log(x) as the natural logarithm ln
        length = CatP.size
        PvIarray = np.empty([length, 3])
        PvIarray[:,0]=np.log(CatI)
        PvIarray[:,1]=np.log(CatP)
        PvIarray[:,2]= CatdP/CatP # Uncertainties in log scale
        # Defining linear function to be fitted by curve_fit
        def lnPvlnI(x, a, b):
            return a*x + b
        # Defining power-law function to be fitted by curve_fit
        def PvItest(x, a, b):
            return b*x**a
        # Fitting power law to data
        if weighted=='true':
            PvIpow, PvIcov = optimize.curve_fit(lnPvlnI, PvIarray[:,0], PvIarray[:,1], sigma=PvIarray[:,2])
            Testpow, Testcov = optimize.curve_fit(PvItest, CatI, CatP, sigma=CatdP)
            # Printing test result
            print('===')
            print('Curve fit test')
            print('Log scale fit: ('+str(PvIpow[0])+','+str(math.exp(PvIpow[1]))+ ')')
            print('Direct power law fit: ('+str(Testpow[0])+','+str(Testpow[1])+ ')')
            print('===')
            print('')
        else:
            PvIpow, PvIcov = optimize.curve_fit(lnPvlnI, PvIarray[:,0], PvIarray[:,1])
            Testpow, Testcov = optimize.curve_fit(PvItest, CatI, CatP)
            # Printing test result
            print('===')
            print('Curve fit test')
            print('Log scale fit: ('+str(PvIpow[0])+','+str(math.exp(PvIpow[1]))+ ')')
            print('Direct power law fit: ('+str(Testpow[0])+','+str(Testpow[1])+ ')')
            print('===')
            print('')
        # Uncertainties on the fit
        PvIerr = np.sqrt(np.diag(PvIcov))
        # Converting the coefficient back to its non-logarithmic value
        PvIcoeff = math.exp(PvIpow[1])
        PvIdcoeff = math.exp(PvIpow[1]) * PvIerr[1]
        # Printing the results of the fit
        print('Power law index: '+str(PvIpow[0])+' ± '+str(PvIerr[0]))
        print('Coefficient: '+str(PvIcoeff)+' ± '+str(PvIdcoeff))
        # Calculating the reduced Chi-Squared value
        if weighted=='false':
            FakeData = lnPvlnI(PvIarray[:,0],PvIpow[0],PvIpow[1])
            Chisq = np.sum((PvIarray[:,0]-FakeData)**2)
            print('Chi-Squared: ' + str(Chisq))
            print('Number of elements: ' + str(PvIarray[:,0].size))
            RedChisq = Chisq/(PvIarray[:,0].size-2) # Reduced Chi-Squared is number of elements minus number of fitted parameters
            print('Reduced Chi-Squared: ' +str(RedChisq))

        if Source != None:
            PvI_plot.text(0.92, 0.85, Source, verticalalignment='bottom', horizontalalignment='right', fontweight='bold')

    # Adding power law fit
    if showfit=='true':
        def PvIfunc(x, A, B):
            return A*x**B
            
        # Creating x values
        PvIfuncX = np.arange(Iscale[0], Iscale[1], 0.1)
        # Calculating the y values
        PvIfuncY = PvIfunc(PvIfuncX,PvIcoeff,PvIpow[0])
        plt.plot(PvIfuncX, PvIfuncY, 'red')
        # Rounding up the fit error
        RoundErrorPvI = round_decimals_up(PvIerr[0],2)
        RoundErrorCoeff = round_decimals_up(PvIdcoeff,1)

        # Adding the fit results as a label
        if weighted=='false':
            Coefficients = '\n'.join((r'$\alpha = {0:.2f} \pm {1:.2f}$'.format(PvIpow[0],RoundErrorPvI),
                                    r'$A = {0:.1f} \pm {1:.1f}$'.format(PvIcoeff,RoundErrorCoeff),
                                    r'$\chi^2_r = {0:.2f}$'.format(RedChisq)))
            PvI_plot.text(0.92, 0.69, Coefficients, verticalalignment='bottom', horizontalalignment='right')
        if weighted=='true':
            Coefficients = '\n'.join((r'$\alpha = {0:.2f} \pm {1:.2f}$'.format(PvIpow[0],RoundErrorPvI),
                                    r'$A = {0:.1f} \pm {1:.1f}$'.format(PvIcoeff,RoundErrorCoeff)))
            PvI_plot.text(0.92, 0.75, Coefficients, verticalalignment='bottom', horizontalalignment='right')

    return PvI_plot

# Function to round error decimals up
#  See https://kodify.net/python/math/round-decimals/
def round_decimals_up(number:float, decimals:int=2):
    """
    Returns a value rounded up to a specific number of decimal places.
    """
    if not isinstance(decimals, int):
        raise TypeError("decimal places must be an integer")
    elif decimals < 0:
        raise ValueError("decimal places has to be 0 or more")
    elif decimals == 0:
        return math.ceil(number)

    factor = 10 ** decimals
    return math.ceil(number * factor) / factor

# =============================================================================
# Polarization fraction as a function of column density
# =============================================================================
def PvN(CatP, CatdP, CatN, Nscale=None, Pscale=None, showfit='false', weighted='true', errorbars='false', FigSize=[5.12,3.84], Source=None):  
    PvN_plot = plt.figure(figsize=(FigSize[0],FigSize[1])) # Creating the figure object
    plt.xlabel('Column Density $N_{H_2}$ (cm$^{-2}$)')
    plt.ylabel('Polarization Fraction $P$ (%)')
    # Plots with and without errorbars 
    if errorbars=='true':
        markers, caps, bars = plt.errorbar(CatN, CatP, yerr=CatdP, fmt ='.', color='black', ecolor='gray', markersize=1, mfc='none',elinewidth=0.9)
        [bar.set_alpha(0.25) for bar in bars]
    else:
        plt.plot(CatN, CatP, '.', color='black', markersize=1, mfc='none')   
    plt.xscale('log')
    axes = plt.gca() # Calling the axes object of the figure
    axes.xaxis.set_ticks_position('both') # Adding ticks to each side
    axes.yaxis.set_ticks_position('both') # Adding ticks to each side
    plt.tight_layout() # Using all available space in the plot window

    # Checking the automated limits for I and P
    xmin, xmax = plt.xlim()
    ymin, ymax = plt.ylim()
    # Setting mininum and maximum values of the I plot
    if Nscale==None:
        Nscale = np.empty(2)
        Nscale[0] = xmin
        Nscale[1] = xmax
    if Pscale==None:
        Pscale = np.empty(2)
        Pscale[0] = ymin
        Pscale[1] = ymax
    # Setting the minimum and maximum values
    plt.xlim(Nscale[0],Nscale[1])
    plt.ylim(Pscale[0],Pscale[1])


    if showfit=='true':
        # Create the array with ln(I), ln(P), and d(ln(P))
        # The numpy package uses log(x) as the natural logarithm ln
        length = CatP.size
        PvNarray = np.empty([length, 3])
        PvNarray[:,0]=np.log(CatN)
        PvNarray[:,1]=np.log(CatP)
        PvNarray[:,2]=CatdP/CatP
        # Defining linear function to be fitted by curve_fit
        def lnPvlnN(x, a, b):
            return a*x + b
        # Fitting power law to data
        if weighted=='true':
            PvNpow, PvNcov = optimize.curve_fit(lnPvlnN, PvNarray[:,0], PvNarray[:,1], sigma=PvNarray[:,2])
        else:
            PvNpow, PvNcov = optimize.curve_fit(lnPvlnN, PvNarray[:,0], PvNarray[:,1])
        # Uncertainties on the fit
        PvNerr = np.sqrt(np.diag(PvNcov))
        # Converting the coefficient back to its non-logarithmic value
        PvNcoeff = math.exp(PvNpow[1])
        PvNdcoeff = math.exp(PvNpow[1]) * PvNerr[1]
        # Printing the results of the fit
        print('Power law index: '+str(PvNpow[0])+' ± '+str(PvNerr[0]))
        print('Coefficient: '+str(PvNcoeff)+' ± '+str(PvNdcoeff))

        if Source != None:
            PvN_plot.text(0.92, 0.85, Source, verticalalignment='bottom', horizontalalignment='right', fontweight='bold')

    # Adding power law fit
    if showfit=='true':
        def PvNfunc(x, A, B):
            return A*x**B
            
        # Creating x values
        PvNfuncX = np.arange(Nscale[0], Nscale[1], 1e20)
        # Calculating the y values
        PvNfuncY = PvNfunc(PvNfuncX,PvNcoeff,PvNpow[0])
        plt.plot(PvNfuncX, PvNfuncY, 'red')
        # Rounding up the fit error
        RoundErrorPvN = round_decimals_up(PvNerr[0],2)
        RoundErrorCoeff = round_decimals_up(PvNdcoeff/1e20,1)
        RoundCoefft = (PvNcoeff/1e20)

        # Adding the fit results as a label
        Coefficients = '\n'.join((r'$\gamma = {0:.2f} \pm {1:.2f}$'.format(PvNpow[0],RoundErrorPvN),
                                 ' '))

        PvN_plot.text(0.92, 0.75, Coefficients, verticalalignment='bottom', horizontalalignment='right')

    return PvN_plot

# =============================================================================
# Polarization fraction as a function of temperature
# =============================================================================
def PvT(CatP, CatT, CatdP=None, Tscale=None, Pscale=None, errorbars='false', FigSize=[5.12,3.84], Source=None, Position='right'):  
    PvT_plot = plt.figure(figsize=(FigSize[0],FigSize[1])) # Creating the figure object
    plt.xlabel('Dust Temperature $T_d$ (K)')
    plt.ylabel('Polarization Fraction $P$ (%)')
    # Plots with and without errorbars 
    if errorbars=='true':
        markers, caps, bars = plt.errorbar(CatT, CatP, yerr=CatdP, fmt ='.', color='black', ecolor='gray', markersize=1, mfc='none',elinewidth=0.9)
        [bar.set_alpha(0.25) for bar in bars]
    else:
        plt.plot(CatT, CatP, '.', color='black', markersize=1, mfc='none')   
    plt.plot(CatT, CatP, '.', color='black', markersize=1, mfc='none')
    #plt.xscale('log')
    axes = plt.gca() # Calling the axes object of the figure
    axes.xaxis.set_ticks_position('both') # Adding ticks to each side
    axes.yaxis.set_ticks_position('both') # Adding ticks to each side
    plt.tight_layout() # Using all available space in the plot window

    # Checking the automated limits for I and P
    xmin, xmax = plt.xlim()
    ymin, ymax = plt.ylim()
    # Setting mininum and maximum values of the I plot
    if Tscale==None:
        Tscale = np.empty(2)
        Tscale[0] = xmin
        Tscale[1] = xmax
    if Pscale==None:
        Pscale = np.empty(2)
        Pscale[0] = ymin
        Pscale[1] = ymax
    # Setting the minimum and maximum values
    plt.xlim(Tscale[0],Tscale[1])
    plt.ylim(Pscale[0],Pscale[1])

    if Source != None:
        if Position == 'right':
            PvT_plot.text(0.92, 0.85, Source, verticalalignment='bottom', horizontalalignment='right', fontweight='bold')
        if Position == 'left':
            PvT_plot.text(0.18, 0.85, Source, verticalalignment='bottom', horizontalalignment='left', fontweight='bold')

    return PvT_plot